/*
 * GMMRecognition.c
 *
 * Created: 31/01/2014 13:31:08
 *  Author: edielson
 */ 

#include "math.h"
#include "GMMRecognition.h"
#include "stdio.h"
#include "stdlib.h"
/*********************************************************
    private constants.
*********************************************************/

//#define PI 2*asin(1)

#define INF 1.5e308
#define PI					3.14159265358979f

/*********************************************************
    private types.
*********************************************************/

typedef struct tag_Mixture
{
	double p;   // mixture weight
	double* m;  // gaussian mean
	double* s; // covariance matrix (only the diagonal elements are stored)
}ttag_Mixture;

typedef struct tag_GMMRecognition
{
	int iNumGaussian;
	int iDimension;
	ttag_Mixture* ptagMixture;
}ttag_GMMRecognition;

/*********************************************************
    private variable.
*********************************************************/

/*********************************************************
    private operations.
*********************************************************/

Bool8	GMMRecognition_loadGMM(Handle hHandle, char* szGMMFilename);
double	GMMRecognition_mixture(double* x, ttag_Mixture* l, int nGaussians, int dim);
double	GMMRecognition_fdp(double*x, double* media, double* var,int ordem);
/*********************************************************
    Operations implementation
*********************************************************/
/*

	Creates the class instance. 

	@param	char* - The GMM file name
		
	@return Handle - The class instance handler 
	
	@author edielsonpf
	@date   12/02/2014
*/

Handle GMMRecognition_new(pInt8 szGMMFilename)
{
	ttag_GMMRecognition* ptagGMMRecog = (ttag_GMMRecognition*)malloc(sizeof(ttag_GMMRecognition));
	if(ptagGMMRecog)
	{
		//if(GMMRecognition_loadGMM((Handle)ptagGMMRecog, szGMMFilename)==False8)
		//{
		//	free(ptagGMMRecog);
		//	ptagGMMRecog = NULL;
		//}		
	}	
	return (Handle)ptagGMMRecog;
}

/*

	Deletes the class instance. 

	@param	Handle - The class instance handler 
		
	@return Bool8 - True8 if succeeded or False8 otherwise 
	
	@author edielsonpf
	@date   12/02/2014
*/

Bool8 GMMRecognition_delete(Handle hHandle)
{
	if(!hHandle) return False8;
	ttag_GMMRecognition* ptagGMMRecog = (ttag_GMMRecognition*)hHandle;
	int iCounterI;
	
	/*Deallocating memory for Mixture*/
	for(iCounterI=0;iCounterI<ptagGMMRecog->iNumGaussian; iCounterI++)
	{
		free(ptagGMMRecog->ptagMixture[iCounterI].m);
		free(ptagGMMRecog->ptagMixture[iCounterI].s);
	}
	free(ptagGMMRecog->ptagMixture);
	return True8;
}

/*

	Loads the GMM model for recognition. 

	@param	Handle - the GMMRecognition class handler
	@param	char* - the GMM file name
			
	@return Bool8 - True8 if succeeded and False8 otherwise 
	
	@author edielsonpf
	@date   12/02/2014
*/
Bool8 GMMRecognition_loadGMM(Handle hHandle, char* szGMMFilename)
{
	if(!szGMMFilename)	return False8;
	if(!hHandle)		return False8;
	
	ttag_GMMRecognition* ptagGMMRecog = (ttag_GMMRecognition*)hHandle;
	Bool8 bResult = True8;
	FILE *ptr;
	int i,j;
	//int Aux;
	
	ptr = fopen(szGMMFilename,"rb");
	fread(&ptagGMMRecog->iNumGaussian,sizeof(int),1,ptr);
	fread(&ptagGMMRecog->iDimension,sizeof(int),1,ptr);
		
	// Allocating memory for GMMs
	ptagGMMRecog->ptagMixture = malloc(sizeof(ttag_Mixture)*(ptagGMMRecog->iNumGaussian));
	for (i=0;i<(ptagGMMRecog->iNumGaussian);i++)
	{
		ptagGMMRecog->ptagMixture[i].m = malloc(sizeof(double)*(ptagGMMRecog->iDimension));
		ptagGMMRecog->ptagMixture[i].s = malloc(sizeof(double)*(ptagGMMRecog->iDimension));
		
		if(ptagGMMRecog->ptagMixture[i].m && ptagGMMRecog->ptagMixture[i].s)
		{
			fread(&ptagGMMRecog->ptagMixture[i].p,sizeof(double),1,ptr);
			for (j=0;j<ptagGMMRecog->iDimension;j++)
				fread(&ptagGMMRecog->ptagMixture[i].m[j],sizeof(double),1,ptr);
			for (j=0;j<ptagGMMRecog->iDimension;j++)
				fread(&ptagGMMRecog->ptagMixture[i].s[j],sizeof(double),1,ptr);	
		}
		else
		{	
			break;
			bResult = False8;
		}
	}
	fclose(ptr);
	return bResult;
}

/*

	Loads the GMM model for recognition. 

	@param	double** pdX- sequence of observation vectors
	@param  int iDim - dimension of each observation vector
	@param  int iNumFrames - number of observations
	@param  ttag_Mixture* ptagMixture - The GMM mixture
	@param  int iNumGaussians - the number of Gaussians in the mixture
				
	@return double - the probability 
	
	@author edielsonpf
	@date   12/02/2014
*/
double GMMRecognition_aposteriori(Handle hHandle, double* pdX, int iDim)
{
	if(!hHandle) return 0;
	ttag_GMMRecognition* ptagGMMRecog = (ttag_GMMRecognition*)hHandle;
	
	double p;		// P(i|X,M)
		
	p = log10(GMMRecognition_mixture(pdX,ptagGMMRecog->ptagMixture,ptagGMMRecog->iNumGaussian,iDim));
	
	return p;
}
/*

	Calculates the P(x[t]|M) probability . 

	@param	double** x- input vector
	@param  ttag_Mixture* l - The Gaussian mixture
	@param  int nGaussian - The number of Gaussians in the mixture
	@param  int dim - the vector dimension
				
	@return double - the probability 
	
	@author edielsonpf
	@date   12/02/2014
*/
double GMMRecognition_mixture(double* x, ttag_Mixture* l, int nGaussians, int dim)
{
	double p; // P(O|M)
	int j;    // counter
	
	p = 0.0;
	for (j=0;j<nGaussians;j++)
		p += l[j].p * GMMRecognition_fdp(x,l[j].m,l[j].s,dim);
	return p;
}
/*

	Calculates the multidimensional Gaussian fdp. 

	@param	double** x- input vector
	@param  double* media - The mean vector
	@param  double* var - the variance vector
	@param  int ordem - the order
				
	@return double - the probability 
	
	@author edielsonpf
	@date   12/02/2014
*/
double GMMRecognition_fdp(double* x, double* media, double* var,int ordem)
{
	// Declaracao das variaveis locais
	double aux=0.0,aux1,aux2 = 0.0; // variaveis auxiliares
	double det; // determinante da matriz de covariancia
	double dif; // var aux p/ o calculo da gaussiana
	int i; // contador
	double prob=0.0; // probabilidade do vetor dada a distribuicao

	// (2*pi)^(ordem/2)
	aux1 = 2.0*PI;
	aux2 = ordem / 2.0;
	aux1 = pow(aux1,aux2);
	//aux1 = exp(aux2*log(aux1));
	
	// Calculando determinante da matriz de covariancia
	det = 1.0;
	for (i=0;i<ordem;i++)
	det *= var[i];

	if (det != 0)
	{
		aux2 = fabs(det);
		aux2 = pow(aux2,0.5);

		for (i=0;i<ordem;i++)
		{
			dif = x[i] - media[i];
			aux += (dif*dif)/var[i];
		}
		aux *= (-0.5);
		aux = exp(aux);

		prob = aux/(aux1*aux2);
	}
	return(prob);
}